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In order to discern the physical nature of many gamma-ray sources in the sky, we must look not only in spectral 
and spatial dimensions, but also understand their temporal variability. However, timing analysis of sources 
with a highly transient nature, such as magnetar bursts, is difficult: standard Fourier techniques developed for 
long-term variability generally observed, for example, from AGN often do not apply. Here, we present newly 
developed timing methods applicable to transient events of all kinds, and show their successful application to 
magnetar bursts observed with Fermi/G-BM. Magnetars are a prime subject for timing studies, thanks to the 
detection of quasi-periodicitics in magnetar Giant Flares and their potential to help shed light on the structure 
of neutron stars. Using state-of-the art statistical techniques, we search for quasi-periodicities (QPOs) in a 
sample of bursts from Soft Gamma Repeater SGR J0501+4516 observed with Fermi/ GBM and provide upper 
limits for potential QPO detections. Additionally, for the first time, we characterise the broadband variability 
behaviour of magnetar bursts and highlight how this new information could provide us with another way to 
probe these mysterious objects. 



1. Introduction 



Neutron stars present the best test cases for ex- 
treme physics in the high-density regime. A long- 
standing problem in neutron star physics is our lack 
of understanding of the neutron star interior, in par- 
ticular, the dense matter equation of state Lattimer 



and Prakash 2007 



The detection of quasi-periodic 
oscillations (QPOs) in the tails of giant flares from 
Soft Gamma Repeaters (SGRs) has opened up the 
possibility of studyi ng neutron s tar interiors using as- 
teroseismology (see Watts|[2012 for a review). 

SGRs exhibit regular bursts in the hard X-rays and 
soft 7-rays (< 100 keV) , and very rare giant flares with 
extremely high isotropic equivalent radiated energy of 
up to 10 46 erg (see e.g. Palmer et al. 2005). Obser- 



vations of persistent soft X-ray counterparts showing 
coherent pulsations with l arge p eriods of 5 — 8 seconds 
Kouveliotou et al. 1998 1999 , and the detection of 



the same periodicities in the tails of the giant flares 
Hurley et al. 1999 Palmer et al.|2005 , suggested that 



SGRs are neutron stars. Their behavior is understood 



within the context of the magnetar model Thompson 



and Duncan 1995 : in this paradigm the SGRs are 
isolated neutron stars with exceptionally strong ex- 
ternal dipole magnetic fields, with internal fields that 



may be as high as 1 10 16 G. Giant flares are pow- 
ered by a catastrophic reordering of the magnetic field 
Woods et al.||2001|. Since this field is coupled to the 



solid crust, Duncan 1998 suggested that such large- 



scale reconfiguration might rupture the crust, creating 
global seismic vibrations that would be visible as pe- 
riodic modulations of the X-ray and 7-ray flux. This 
idea was confirmed by the detection of QPOs in the 
expected range of frequencies (~ 10 — 1000 Hz) in the 
tails of giant flares from two different magnetars Is- 
rael et al.||2"005] |Strohmayer and Watts] [20051 [2006] ' 



Watts and Strohmayer 2006 . If the QPO frcquen 



cies can be reliably identified with particular global 
seismic modes of the neutron star, then they can in 
principle be used to constrain both the equation of 
state and the interior magnetic field. 

A major obstruction to this field of research is the 
sparsity of data. Since the launch of the first X- 
ray and 7-ray instruments, only three giant flares 



1 Supported by period and pe- 

riod derivative measurements; see 

http : //www. physics. mcgill.ca/~pulsar/magnetar/main. html 
for an up-to-date reference list on magnetar spin-down 
properties 
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have been observed, with just two having data with 
a sufficient time resolution to detect QPOs. In try- 
ing to overcome this lack of observational constraints, 
it is therefore a reasonable approach to turn to the 
much more numerous short SGR bursts with lengths 
of usually less than a second and luminosities around 
10 40 ergs -1 . Hundreds of SGR bursts have now been 
observed from many magnetars 2 . 

To date there has been no systematic search for pe- 
riodic features in the lightcurves of the SGR bursts. 
A search for QPOs in a period of enhanced emission 
with multiple bursts (a 'burst storm'), from the mag- 
netar SGR J1550-5418, carried out using data from 
the Fermi Gamma-ray Burst Monitor (GBM), found 
no significant signals |Kaneko et al.|2010 



and Ibrahim 2010 



El-Mezeini 



Rossi X-ray 



searched a subset of 
Timing Explorer data from SGR 1806-20 for periodic 
features and found some tentative signals: however, 
there are several points of concern with regard to their 
methodology explained in full in the Appendix of |Hup-| 



penkothen et al. 2013 . 



Searching for QPOs in transient light curves is a 
non-trivial task. Standard methods involve Fourier 
analysis, more specifically the periodogram, defined 
as the squared amplitude of the Fourier transform of 
the light curve. The periodogram has several advan- 
tageous statistical properties: Poisson noise prevalent 
in light curves from photon counting experiments re- 
sults in a flat periodogram with a well-known statis- 
tical distribution, making the detection of periodic 
and quasi-periodic features (as outliers above that 
distribution) tractable. However, the detectability of 
QPOs changes significantly in the presence of corre- 
lated noise processes and the transient properties of 
the bursts we are concerned with here. The very na- 
ture of a transient event - it has a start, one or more 
peaks, and an end - complicates the analysis procedure 
and introduces additional sources of uncertainty, es- 
pecially in the low-frequency part of the periodogram. 
For transient events where the shape of the burst en- 
velope is known, many problems arising from the non- 
stationarity can be solved either analytically Guidorzi 
2011 or via Monte Carlo simulations Fox et al.|2001 



However, many astrophysical transients such as mag- 
netar bursts do not show a well-behaved burst light 
curve that is easily reproducible by a simple function. 
This in itself can be interesting, aside from searching 
for QPOs: the different burst envelope shapes must be 
created by a physical process in the source, either in 
the form of noise processes or non-stochastic emission 
processes, and characterising the differences may tell 



20061 



e-g- 

|mT 



:ghctti 



Woods 
[2008T 



for 



and Thompson| 



overviews 

http:/ /f64. nsstc.nasa.gov/gbm/science/magnetars/ for 



us more about the emission processes at work. The 
methods and analysis presented in the following sum- 
marise more extensive work laid out in |Huppenkothen| 



et al. 2013 . We refer the reader to that paper for 



more details. 



2. Variability Analysis 

2.1. Monte Carlo Simulations of Light 
Curves - Advantages and Shortcomings 

Monte Carlo simulations of light curves are a stan- 
dard tool in timing analysis (see for example Fox et al. 
2001). The underlying idea is simple: one fits an em- 



collection of SGR bursts observed with Fermi/ GBM 



pirically derived (or physically motivated) function to 
the burst light curve. One then generates a large 
number of realizations of that burst profile, includ- 
ing appropriate sources of noise, such as Poisson pho- 
ton counting noise. The periodograms computed from 
these fake light curves form a basis against which to 
compare the periodogram of the real data. For each 
frequency bin, a distribution of powers is produced. 
Comparing the observed power in each bin with the 
distribution of simulated powers in the same bin al- 
lows us to make a statement about the probability of 
the observed power in a particular bin being due to 
a noise process: if the observed power in a particular 
bin is a high outlier compared with the distribution 
of simulated powers in that bin, then the probability 
of observing the data under the (null) hypothesis of a 
noise process is 1/N, where N is the number of simu- 
lations performed. If N is large, the observed outlier 
is unlikely to be produced by the noise process alone. 

The Monte Carlo method outlined above is versa- 
tile and powerful, but it has limitations. The most 
important limitation comes from our lack of knowl- 
edge of the underpinning physical processes producing 
the observed light curve. Only if the null hypothesis 
accurately reflects the data - apart from the (quasi- 
periodic signal for which we would like to test - is 
the test meaningful. If important effects that distort 
either shape or distribution of the powers are missed, 
then the predictions made will not be accurate, lead- 
ing to either spurious detections or real signals not 
being found. 

More often than not, especially in the case of short 
magnetar bursts, we do not have complete informa- 
tion about the emission mechanism. Short magnetar 
bursts are extremely diverse, varying in light curve 
shape as well as burst intensi ty a nd duration (see, fo r 
example, G6gu§ et al. |1999 and Gogiis, et al. 2000 ) . 
There is a fundamental degeneracy in the problem: 
which features to we consider to be part of the burst 
envelope, and which a (potential) feature which we 
do not include in the light curve template? Stochas- 
tic processes correlated on different time scales can 
mimic a quasi-periodic process to the human eye, but 
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clearly, this is not what we wish to detect. Not taking 
the presence of these features into account can lead 
to a large number of false positive detections. On the 
other hand, overcorrecting for features may lead us 
to detect no features at all. Without detailed knowl- 
edge of the burst emission processes, it is impossible to 
build a reliable model for individual burst light curves 
to test agains when searching for QPOs. We therefore 
advocate a different method, based on an empirical 
model of the periodogram. 



2.2. Modeling the Periodogram 

We give a very short overview of the general prin- 
ciples our method employs. Details can be found in 



Huppenkothen et al. 2013] and Vaughan 2010 The 



Bayesian procedure we employ in modeling the peri- 
odogram has three parts: (a) find the preferred broad- 
band noise model to represent the low-frequency part 
of the periodogram, (b) search the periodogram for 
the highest outlier and compare this outlier to those 
distributed by pure broadband noise to find narrow 
features, (c) search for QPOs in the data, using binned 
data as well as an identical approach for the model 
selection in the first step. A step-by-step description 
can be found in the Appendix of Huppenkothen et al. 



2013 



Every step in the analysis follows the same logic: 
assume a null hypothesis and an alternative hypothe- 
sis, compute statistics to summarise the data-model 
fits for the two different models, generate a sam- 
ple from this null hypothesis using a Markov Chain 
Monte Carlo (MCMC) approach, more specifically the 
MCMC code emcee [Foreman-Mackey et "aL] |2012 



One can then compare the distribution of the relevant 
statistic derived from the sample generated from the 
null hypothesis to the observed value of that statistic. 
If the observed value lies in the high-end tail of the 
distribution, then it is an outlier with respect to the 
null hypothesis. 

Since the entire procedure rests on the correct 
choice of broadband model, this is the first step of 
the analysis. The data are fitted with two continuum 
noise models, which, by definition of the likelihood 
ratio test, are required to be nested. The likelihood 
ratio is the statistic we use to decide which model 
is preferred by the data. We simulate a large num- 
ber of fake periodograms from parameter sets drawn 
from the posterior distribution of parameters, as ap- 
proximated by a large number of MCMC simulations. 
Then these fake periodograms are fit with both mod- 
els again to build a distribution of likelihood ratios 
from the simple model. We can compute the tail-area 
probability (p-value) of the observed likelihood ratio 
to be typical of the distribution (equivalent to asking 
whether the observed data is sufficiently described by 
the simpler model) by integrating over the tail of the 



distribution. If this probability is lower than a cho- 
sen significance threshold, then the data is more likely 
to be drawn from the more complex model hypothe- 
sis, which should then be adopted for the rest of the 
analysis. 

We extensively tested our method on synthetic 
data generated with known parameters. These sim- 
ulations confirm that in the limit of white noise, 
our method matches the predictions from standard 
Fourier analysis. Furthermore, for more complicated 
light curves, involving red noise and a burst envelope, 
our method is conservative in nature at low frequen- 
cies, where burst envelope and red noise dominate, 
but approaches the white noise predictions for high 
frequencie s. The full description of th ese simulations 
is given in 



Huppenkothen et al. 2013 



3. Observations 

To test our methodology, we applied it to a small 
sample of bursts from SCR J0501+4516. Fermi/GSM 
triggered 26 times on this source between 2008 August 
22 and 2008 September 03, observing 29 bursts. Two 
of these (burst IDs 080824054 and 080825200) had 
saturated parts, and were therefore excluded from the 
analysis due to the rather complicated effects satura- 
tion can have on periodograms. Following |Lin et aL] 



2011 



we used only Nal detectors with a source an- 
gle smaller than 50deg for each of the 24 triggered 
and 3 untriggered bursts. We use high-time resolution 
time-tagged event (TTE) data, for which time and 
energy for each individual photon is recorded. The 
data were barycentered and channels converted to the 
mid-energy of each energy bin. The observations were 
then energy-selected to include only counts between 
8 and 100 keV. The lower limit to the energy is set 
by the detector response Meegan et al.|2009 , the up- 
per limit was found by inspecting energy-resolved light 
curves and finding no source counts above 100 keV 
(the counts were consistent with the Poisson distribu- 
tion expected from counting noise). Burst start times 



and lengths (T90 durations) were taken from Lin et al. 
[2011] , and are summarized in Table 1 of that paper. 
We added 20% of the burst duration to both ends of 
the burst in order to ensure that we caught the entire 
burst. A selection of six bursts is shown in Figure [TJ 
to emphasize the diversity of burst morphologies we 
encounter. 



4. Results 

We computed light curves and periodograms for 
all 27 bursts. In each case, we produced a light 
curve by binning the TTE data to a time resolu- 
tion of l/2v Nyquist = 1.22 x 10~ 4 s, corresponding 
to a Nyquist frequency of t'Nyquist = 4096 Hz. We 
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Time since trigger [s] 



Figure 1: Light curves of six example bursts from the magnetar SGR J0501 + 4516 recorded by Fermi/GBM. We 
combined data from all Nal detectors with source angles smaller than 50 degrees to the source. The time resolution 
corresponds to 0.005 seconds. Note the strong component of aperiodic variability after the main burst in 080823478 and 
the differences in peak count rate by almost one order of magnitude between the upper three bursts and the lower three. 




Time since trigger [s] log(Frequency) [Hz] 

Figure 2: Fermi/ GBM observation of burst 080823847a from SGR J0501 + 4516. Left: light curve with a time 
resolution of 0.002 seconds. Structure in the burst profile is clearly visible. Right: unbinned periodogram (blue) and 
two examples of binned (magenta: 16 Hz binning; orange: 65 Hz binning) periodograms for this burst. There is a 
feature in the periodogram around 30 Hz (leftmost arrow), which is by itself not significant. However, significant 
features (arrows 2-5) are all at integer multiples of this frequency (within the uncertainty imposed by the frequency 
resolution), indicating the presence of harmonics at 150 Hz, 300 Hz, 900 Hz and 2100 Hz. 



chose the time resolution based on the Nyquist fre- 
quency of interest: we do not expect any signals above 



4000 Hz from neutron star seismic oscillations McDer- 
mott et ai"||1988| . We search both the unbinned peri- 



odogram as well as the same periodogram binned to 
integer multiples (3, 5, 7, 10, 15, 20, 30, 70, 100, 200, 
300, 500 and 700) of the frequency resolution of that 
burst. 



4.1 . Search for QPOs 

None of the 27 bursts shows periodicities of any 
significance in any of the unbinned periodograms. 
The highest data/model outlier significance is seen in 
burst bn080823847a (see Figure [| for a light curve 
and periodogram), with a posterior p- value £>(Tr) = 
0.11 ± 0.01, at frequency ^ max = 4057 Hz, well below 
the power required to reach the detection threshold 
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corresponding to a posterior p- value of 5%. However, 
the same burst shows significant signals from ~300 Hz 
to ~2100 Hz in the several binned periodograms (bin- 
ning frequencies between 16 Hz and 1583 Hz). The 
most significant signals are for binning frequencies 95 
Hz and 158 Hz at a frequency of ~2100 Hz with a 
posterior p- value < 2.0 x 10~ 5 . Note that p- values 
quoted there are corrected for the number of frequen- 
cies searched, but neither for the number of bursts 
searched nor the number of binned spectra searched 
for each burst. While the former is straightforward (a 
simple multiplication factor of 27 for the number of 
bursts searched), the latter is more complicated, ow- 
ing to the fact that searching different binnings for a 
single periodogram does not result in independent tri- 
als. The most conservative assumption is to consider 
them independent, including another multiplication 
factor equal to the number of binnings searched (here: 
9). This would rule out all but the two signals with 
frequency bins of 95 and 158 Hz, which remain signifi- 
cant even after a correction for the number of trials. A 
possible interpretation of the observed signals follows 
in Section [5] 



4.2. Broadband Variability 

The broadband variability observed in the bursts is 
not just a nuisance when searching for (quasi-) peri- 
odicities, but is of interest in its own right: it shows 
that something is varying in the source, although not 
periodically. Almost all bursts in the sample are well- 
modeled by a simple power law. The distribution of 
indices ranges from 1.7 to 4.3 and peaks around 2.5, 
which is higher than commonly seen for example in 



Ga mma-ray bursts (see e .g. Beloborodov et al. 2000 
and Guidorzi et al.|[2012| . Only two bursts required 
the more complex model. While these were not the 
longest bursts, they had the highest fluence (except for 
the excluded saturated bursts) , indicating a potential 
correlation between power spectral shape and burst 
fluence. A reliable characterisation of the broadband 
properties will be deferred to a future paper involving 
a larger sample of bursts. 



5. Discussion 

Magnetar bursts are a potential window into the in- 
terior of neutron stars, via the oscillations measured 
in magnetar giant flares. Finding analogous signals 
in the wealth of short SGR bursts, however, poses 
something of a challenge. We have shown that tim- 
ing analysis of astrophysical transients is a non-trivial 
problem. Standard Fourier techniques are insufficient 
for phenomena with diverse light curve morphologies, 
especially when involving correlated noise processes. 



Monte Carlo simulations of light curves fail to be pre- 
dictive when there is no precise knowledge of the un- 
derlying burst light curve: there is a degeneracy be- 
tween the overall, aperiodic burst shape, a potential 
red noise component, and the very thing we would like 
to measure: a QPO. 

In the absence of better knowledge about the emis- 
sion processes in magnetar bursts, we advocate a con- 
servative Bayesian method that models the burst light 
curve as a pure red noise process, at the cost that 
weak signals are likely missed. This is the greatest 
weakness of our approach. Even strong signals may 
be undetectable at low frequencies, where burst en- 
velope and red noise dominate. This limitation is in 
part not only due to restrictions of our method, but 
also to the short lengths of the SGR bursts, where 
at these frequencies only one or two cycles may be 
seen in the light curve. However, at frequencies close 
to and above 100 Hz, sensitivities approach the white 
noise limit, which is strongly dependent on the num- 
ber of photons from a particular burst. Thus, for a 
bright burst with good count statistics, sensitivities 
are quite constraining, down to less than 10 %. This 
is comparable to what was observed in giant flares: for 
example, a QPO at 93 Hz, as seen in the 2004 flare, at 
roughly 10 % rms am plitude Israel et al.|2005 Watts 



and Strohmayer 2006 , should be detectable in at least 



the brightest bursts of our sample. 

However, QPOs in SGR bursts may be less strong 
than in the giant flares, owing to the lower energy 
injected in SGR bursts, and hence more likely to be 
misclassified as non-detections, if their fractional rms 
amplitudes fall below 5 %. 

The burst 080823847a presents an interesting case 
that illustrates the limits of a pure signal-processing 
approach to the timing analysis shown here. The na- 
ture of the significant signals is at present unclear. 
They are possibly harmonics of a lower-frequency sig- 
nal around 30 Hz, corresponding to a timescale of 
r = l/v = 33 ms. This timescale roughly corresponds 
to the two sharp peaks seen in the burst light curve in 
Figure [2] (left side). Whether we consider this to be a 
QPO atop a burst envelope or not cannot be answered 
from Fourier analysis alone; it becomes a matter of 
interpretation and prior knowledge. At present, with- 
out any knowledge about emission processes and the 
kind of light curve they produce, it is impossible to 
distinguish whether the two sharp peaks are indeed a 
heavily damped QPO, or simply a chance occurrence 
of red noise features, thus we choose the conservative 
approach and interpret the observed feature as part 
of a noise process. 

We wish to note that methods presented here, while 
developed with SGR bursts in mind, are by no means 
limited to magnetars. They are applicable in fairly 
general circumstances, for any light curve that is 
phenomenologically similar to what we observe from 
magnetars: highly variable, transient events with 
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complex light curves. This includes, for example, 
other known transients such as gamma-ray bursts 
(GRBs), tidal disruption events and supernova light 
curves. 
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